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ABSTRACT 

We present a perturbation theory for studying the instabihties of non-axisymmetric 
gaseous discs. We perturb the dynamical equations of self-gravitating fluids in the 
vicinity of a non-axisymmetric equilibrium, and expand the perturbed physical quan- 
tities in terms of a complete basis set and a small non-axisymmetry parameter e. 
We then derive a linear eigenvalue problem in matrix form, and determine the pat- 
tern speed, growth rate and mode shapes of the first three unstable modes. In non- 
axisymmetric discs the amplitude and phase angle of travelling waves are functions of 
both the radius R and the azimuthal angle (f>. This is due to the interaction of different 
wave components in the response spectrum. We demonstrate that wave interaction in 
unstable discs with small initial asymmetries, can develop dense clumps during the 
phase of exponential growth. Local clumps, which occur on the major spiral arms, can 
constitute seeds of gas giant planets in accretion discs. 

Key vifords: accretion discs - hydrodynamics - instabilities - galaxies: structure - 
galaxies: spiral 



1 INTRODUCTION 

Unstable density waves inside self-gravitating gaseous media 
play an important role in structuring of spiral galaxies and 
accretion discs. Prior works in the literature deal with the 
instability of differentially rotating discs, which are initially 
circular. Aoki et al. (1979) carried out a global analysis of 
spiral modes in the soft-centred, gaseous Kuzmin (1956) disc 
using Hunter's (1963) device. Lemos et al. (1991) studied the 
unstable axisymmetric modes of the more puzzling scale- free 
discs and solved the dispersion relation to determine the suf- 
ficient conditions for ring formation and for the existence of 
breathing modes. Recently Goodman & Evans (1999, here- 
after GE) constructed the normal modes of the scale-free 
Mestel (1963) disc in the Cowling approximation. For self- 
gravitating modes they computed the ratio of growth rate to 
pattern speed by solving a recurrence relation in the Mellin- 
transform space. Density waves in the presence of an exter- 
nal potential field, like that of a massive central object, have 
also been extensively studied using semi-analytical and finite 
difference techniques (e.g., Papaloizou & Lin 1989; Savonije 
& Heemskerk 1990; Papaloizou & Savonije 1991). 

The assumption of axisymmetry of the equilibrium 
(or initial) state, however, imposes significant simplifica- 
tions on the solution of dynamical equations. In fact, un- 
stable modes of circular discs are decoupled because the 



* naser_ma@iasbs.ac.ir 
t mjalali@sharif.edu 



only 0-dependent term, exp{im(j}) (i = V— 1), is factored 
out of the linearized Poisson and hydrodynamic equations. 
Here (j) is the azimuthal angle and m is the wave number. 
Such a factorization is impossible if the unperturbed state 
is non-axisymmetric. Azimuthal modes of perturbed non- 
axisymmetric discs can interact with the Fourier compo- 
nents of the equilibrium state and result in more complex 
patterns. 

In this paper we attempt to formulate the instability 
problem of non-axisymmetric discs based on a perturbation 
theory, and present a mechanism for the generation of dense 
clumps through wave interaction in linear regime. That oc- 
curs when a non-axisymmetric disc is destabilized. We apply 
our formulation to the simple discs of Jalali & Abolghasemi 
(2002, hereafter JA), which reduce to Mestel's disc in ax- 
isymmetric limit. We have chosen JA models as our unper- 
turbed state because of two reasons. Firstly, closed form 
(exact) expressions are available for physical quantities like 
the potential-density pair, pressure, sound speed and veloc- 
ity components. Secondly, both the potential and the surface 
density have a single Fourier component. This minimizes the 
number of interacting terms in our subsequent analysis. 

We review steady-state, non-axisymmetric, scale-free 
discs of JA in section 2 and introduce a cutout func- 
tion that helps us handle the singularity at the cen- 
tre. We present a first-order perturbation theory for non- 
axisymmetric gaseous discs in section 3, where we also in- 
troduce a basis set by which we expand the perturbed phys- 
ical quantities. We then derive an eigenvalue problem for 
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Figure 1. The region under the sohd hue is tlie allowable zone 
in the parameter space of a bisymmetric disc with n = 2. 



determining the normal modes of non-axisymmetric discs 
and discuss on the algorithms used for solving the result- 
ing eigensystem. In sections 4 and 5, we study the unstable 
modes of circular and non-axisymmetric discs for the fun- 
damental wave numbers m = 0, 1 and 2, and discuss on the 
results and their physical implications in section 6. 



2 NON-AXISYMMETRIC GASEOUS DISCS IN 
EQUILIBRIUM 

There are few papers in the literature that address non- 
axisymmetric equilibria of gaseous discs (Syer & Tremaine 
1996; Galli ct al. 2001; JA). In this section wc briefly review 
one of exact non-ax;isymmetric models of JA whose surface 
density consists of a simple harmonic function of the az- 
imuthal angle. 

2.1 Basic derivations 

We define {R, (f) as the usual polar coordinates and assume a 
sclf-gravitating, inviscid, compressible fluid in a steady two 
dimensional state so that the surface density distribution 
and its self-consistent potential are scale-free as 
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Here G is the constant of gravitation, and Ro and Es are 
the normalizing constants of the length and surface density, 
respectively. Our discs arc non-axisymmetric because E and 
$ depend on the azimuthal angle 0. Due to the continuity of 
physical quantities in the azimuthal direction, the functions 
f{4>) and g{4>) are 27r-periodic in 0. For the sake of simplicity, 
we work with non-axisymmetric discs of the form 



/ {4>) = 1 + € cos n4 



(3) 



where < e < 1. The 0-dependent part of the potential 
reads 

2 

g{<j>) = — ecosncf). (4) 

For e = 0, our discs reduce to the isothermal, gaseous Mes- 
tel disc. Let us respectively define co, Vc and M=Vc/co as 
the sound speed, velocity of circular orbits and Mach num- 
ber in Mestel's disc. The components of the velocity field 



associated with (1) and (2) are obtained after solving the 
hydrodynamic equations as (JA) 
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The sound speed c{4>) and the pressure distribution P{4>) 

are then given by 



c{(j>) = Co 
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c reduces to co in the limit of e = 0. From (2) and (8), and 
equation (2) in JA, one can verify that Vq—Vc + Cg. Sound 
waves are determinate if the pressure is positive, i.e., > 0. 
A necessary and sufficient condition for the positiveness of 
P is 
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Here Amin is the niininmm of A{4>) over the interval < (?!) < 
27r. Inequality (10) restricts model parameters. JA models 
allow for transonic flows because the sound speed depends 
on (/) and it varies along a given streamline. In principle, 
transition from subsonic to supersonic speeds can take place 
smoothly but the reverse phenomenon is impossible without 
a shock wave. To avoid shock waves, we follow JA and ex- 
clude transonic flows from our calculations and investigate 
supersonic and subsonic discs. 

In this paper we choose a bisymmetric equilibrium with 
n = 2. The allowable region of the parameter space (M, e) 
is determined based on two requirements. Firstly the square 
of the sound speed must be positive. Secondly, transonic 
flows arc excluded. Fig. 1 shows the admissible zone of the 
parameter space. Note that for e > the Mach number is a 
function of the azimuthal angle: 



M(M,e, 



(11) 



In the limit of circular discs (e = 0) we obtain M = M as 
expected. 

2.2 Active surface density 

We are studying scale-free gaseous discs whose surface den- 
sity and potential diverge near the center according to (1) 
and (2). The linear stability problem of scale-free discs is 
ill-posed unless one decides on the way that inward travel- 
ing waves are reflected off the cusp. Several methods have 
been devised to tackle this problem. For scale-free stellar 
discs Zang (1976) (sec also Evans & Read 1998a, b) used 
an inner cutout, which carves out the density function and 
deactivates the central region of the stellar system. In this 
way, stars with high orbital frequencies do not participate in 
perturbations. For the gaseous Mcstcl disc GE used a hard 
internal boundary that fixes the phase with which waves are 
reflected from the centre. In this paper we apply Evans & 
Read's (1998a) method to our cuspy gaseous discs of §2.1 
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and deactivate some parts of the disc mass by the cutout 
function 



H{R) 



Rr, 



We then analyze the perturbations of the active density 



(12) 



Eact(i?) =S(i?)//(i?). 



(13) 



Here Rin = aRo and Tiout = Ro/^ are the inner and outer 
cutout radii, respectively. Since the cutout is not applied 
to the potential, the immobile mass will resemble a rigid 
bulge/halo component that does not participate in the disc 
dynamics, but it contributes to the force field. Instability 
properties are not sensitive to the choice of Rout ■ We apply 
the outer cutout for computational convenience and use a 
largo -Rout far beyond the region that instabilities occur. The 
nature of functions by which we will expand perturbed fields, 
constrains the feasible values of iVin and Nout ■ We have used 
iVin= iVout = 2 throughout this paper. 



3 PERTURBATION THEORY OF 
NON-AXISYMMETRIC DISCS 

We define u = (wi,W2,M3)^ = {T,,vr,v^)^ and consider the 
first order Eulerian perturbations of physical quantities as 

u — »• uo{R,(j)) + ui{R,(l),t), 

$ ^ MRA) + MRA,t), (14) 



c^(0)So(J^, 0) + c\(t>)T.i{R, ct>, t). 



The superscript T denotes transpose and the subscripts 
and 1 refer to the equilibrium and perturbed states, 
respectively. We assume that the perturbed quantities, 
5i s ($1 , El , vm , v^i)^ , are composed of a discrete set of nor- 
mal modes, each of the form 



— H" iSjTj,, i — 1, 



(15) 
(16) 



where uJm, and s-m are respectively the cigcnfrequency, 
pattern speed and growth rate of the mth mode. It is re- 
marked that the real part of (5i,m gives the physical solution. 

We use the equilibrium fields of §2 and substitute from 
(14) into the continuity and momentum equations of self- 
gravitating fluids [equations (lE-3) and (6-20) in Binney 
& Tremaine 1987]. After collecting the first-order terms in 
perturbed variables, we obtain 



L • ui 



(17) 



with = dl'd[i. The elements of the linear operator L = 
[Lij], (j, j = 1, 2, 3), have been given in Appendix A. $i and 
El should also satisfy Poisson's integral 
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Modes of non-axisymmetric discs bifurcate from those 
of circular discs once we let e grow from zero. Therefore, we 
seek for a perturbation solution for 5i,m. in terms of the non- 
axisymmetry parameter e. Moreover, continuity of physical 
quantities in the azimuthal direction implies the periodicity 



conditions 

uoiR,(l> + 2n) = uo(J?,0), 
S["^\R,<P + 2n,e) = <5^'(7?,0,e) 



(19) 
(20) 



These requirements suggest to represent uo by the Fourier 
series 



Uo(J?,0) = 



J2 ^'"'e'" 



'uo,p{R), 



and propose the Ansatz 



(21) 



(22) 



Using equations (21) and (22) and substituting from (15) in 
(17), we obtain the determining equation of the mth mode 
as 



l—np) 



R 

Co 



^(Q,dR¥^+"'\i{m + nl)R-^$'r^"'^) , (23) 
Z = 0,±l,±2,..., 



which is a system of linear ordinary differential equations 
with respect to R. Here DpS are linear differential operators 
(they are square matrices of dimension 3x3) whose ele- 
ments have been given in Appendix B for < |p| < 3. The 
second term on the left side of (23) shows that how different 
azimuthal waves interact. The interacting terms rapidly de- 
cay proportional to e'*'' for e <C 1. In the axisymmetric limit 
(e = 0), equation (15) reduces to 



5i,m{RA,t) = e-""-*e'""^d["'\R), 



(24) 



so that LOm is associated with a single azimuthal wave com- 
ponent. Modal decomposition of this kind is impossible for 
initially non-ax:isymmetric discs, for the equilibrium fields 
include the powers of e*'"'* that cause wave interaction. 



3.1 Choice of basis functions 

One way for solving (18) and (23) is through expanding 
S[''\r) [see equation (22)] in terms of some basis, preferably 
biorthogonal, functions. Glutton-Brock (1972) found such a 
set, which was then recalculated by Aoki & lye (1978) in 
terms of associated Legendre functions. We follow them and 
expand the functions 5^*' (R) as 



j=\k\ 



u[''\r) = ^pl'=l(OW.X, 



(0, (25) 



(26) 



where X=(a*', 6*, cf)^ is a to-be-determined vector of com- 
plex coefficients and 

W = diag[E,wi(0,icoW2(0,coW2(0], (27) 

'i-^v ... _ M-r 



- Rc 
^ = WTr^' 
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(5) are the normalized associated Legendre functions 
defined by 



0--|fc|)! (2j + l) 

2{j + \k\y. 



(28) 



Wo set the scale length of Clutton-Brock functions, Rc, to 
the normalizing length of (1) and (2). i.e., Rc = -Ro- The 
criterion in choosing the weight function U'2(C) was the con- 
vergence of proceeding definite integrals that constitute the 
coefficients of the matrix A in (29). If fc 7^ 0, the functions 
Pj*'(^) give the factor [(1 +C)/2]I''I/^ Therefore, the expan- 
sions used for velocity components are not singular at the 
origin (when ^ — » —1). For fe = 0, the perturbed velocity 
components become singular for ^ = — 1. Nonetheless, the 
energy and angular momentum are not singular there and 
the continuity equation is satisfied. 



3.2 Matrix formulation 

We now substitute from (25) and (26) into (23), left-multiply 
both sides of the resulting equation by — iPj'^'w^^ (fc = 
0, 1, 2, • • •), and integrate over < ii < 00. Using the or- 
thogonality of the associated Legendre functions, we obtain 
a system of linear equations for the coefficients aj, 6* and 



A.x = 5^x, 

Co 



(29) 



X — (■ ■ ■ ) 3-m-t-n: bm-t-n? ^m-\-n: ^m; bm? Cm, ^m—n^ • • •) ' 

a-r = {a\r\,a'\r\+i,a'\r\+2, ■■■) J 

hr = (6[r|,^[r|-|-li^[T-H-2j---) I 
C|r|) + C|^|_|_2, ■ ■ ■) . 

It is seen that uJmRo/co is an eigenvalue of A and X is its 
corresponding eigenvector. We find that the elements of A 
are linear combinations of the following definite integrals 
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where r and s are integer numbers of the form m + np' and 
m -|- np", respectively. The integrals in (30) are evaluated 
using Gaussian quadratures. We transform A to Hessenberg 
form and obtain its corresponding QR transform (Press et 



al. 1997). We then compute the eigenvalues and eigenvectors 
of the linear system (29). Consequently, the surface density 
of the mth mode is determined from (22) and (26) as 



Ei,m(il,<?i,f,e) = e'™k™(i?,<^,f) 



X cos[m4>-Q,mt + Q,n{R, <f), e)]. 



Arr^iR, 0, e) = V E*'"' {R, <p, e)Er ' {R, 0, e), 



(m) / 



Qm{R, (j>, e) = arctan 
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(m) 



(31) 
(32) 

(33) 



where a bar denotes complex conjugate. One can expand 
Am and 6m in the Taylor series of e as 



Y^(R,(l>,t) 



ym,o(i?)+^e'Fw(i?, 



r = (^,e),(34) 



so that Ym.i are 27r-periodic in 4>. The leading zcroth order 
terms of (34), Ym,o{R), determine the overall shape of the 
perturbed density, which is usually an m-fold trailing spiral. 
The higher-order terms (Z > 0) lead to density fluctuations 
along the major spiral arm(s). 

Both A and X arc infinite dimensional. However, the 
non-axisymmetry parameter e is small and its powers, which 
serve as the coefficients of interacting terms, fall off rapidly. 
Therefore, it is reasonable to truncate the series of the inter- 
acting modes in (23) at some \l — p\ = L. We also truncate 
the series in (25) and (26) at jmax = J ■ We begin with L = 1 
and J = 20, and increase them until we acquire a relative 
accuracy of 1% in computing oJm- 

In the axisymmetric limit, A becomes a block diagonal 
matrix and we obtain the following eigenvalue problem for 
the mth mode 



(35) 



,c7), 



■D/ ^ rr <^mRo „ 
B(m) • Z = Z, 

CO 

= [ai ,...,aj ,bi 

with B(m) being a single 3 J x 3J block of A. The perturbed 

surface density corresponding to (35) becomes 

Ei,m(7?, 0, t) = e""'*^™,o(ii) cos [mct)-nmt+em,o{R)] ■ (36) 

The matrix B has 3J eigenvalues, which lie on several 
branches in the complex (n^,, Sm)-plane. However, eigen- 
values of most branches diverge as J is increased. They are 
virtual roots of the determinantal equation 

Ro 



B(m) 



I 



Co 



0, 



(37) 



with I being the unit matrix. Physical eigenvalues of un- 
stable discs (if they exist) belong to a convergent complex 
branch. The convergent eigenvalue with the largest growth 
rate corresponds to the fundamental eigenmode. 

The ideal scale-free limit will be attained if we let both 
a and 7 tend to zero. Our calculations show that the sta- 
bility/instability properties do not change substantially by 
decreasing 7. The choice of 7 = 0.05 gives accurate and ro- 
bust results in most cases. So we need to play with a. Since 
the scale length of Clutton-Brock functions is fixed, by de- 
creasing a we require too many terms in the expansion of 
5['''' (R) to get the series converged. The smallest value of a 
that we could approach, while assuring computational accu- 
racy, was a = 0.03. 
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Figure 2. Eigenvalue spectra of circular discs for a = 0.5. Squares and filled circles correspond to J=50 and .7=55, respectively. The 
converged fundamental eigenfrequency has been shown by uic- Top, middle and bottom rows correspond to m = 0, 1 and 2, respectively. 
Left, middle and right columns are associated with M=0.7, 2.5 and 5, respectively. 
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Figure 3. Variation of the fundamental eigenvalue of circular discs versus M. Solid, dashed and dotted lines correspond to a=0.03, 0.1 
and 0.5, respectively, (a) The growth rate so of the m = mode. For this mode we have f2o = 0. (b) The variation of arg(c<Ji) versus M. 
The argument of iO\ abruptly drops to zero at a critical M. This critical value tends to 1 as a is decreased, (c) Same as Panel 6 but for 
m = 2. The argument of L02 monotonically tends to zero by decreasing M. 



6 Naser M. Asghari and Mir Abbas Jalali 





Figure 4. Top row. unstable modes of a circular disc with a = 0.5 and M = 5. The corresponding absolute values of aj' (m = 0, 1, 2) 
have been plotted versus j below each mode shape. The corotation radius is iicR = 0.972 and iicR = 0.943 for m = 1 and 2, respectively. 



4 MODES OF CIRCULAR DISCS 

To this end, we study instabilities of cutout Mestel discs, 
which are obtained from the models of §2 by setting e = 0. 
Calculation of unstable modes of circular discs is straight- 
forward because the equilibrium fields do not depend on the 
azimuthal angle 4>. We set m = 0, 1 and 2 in (35) and cal- 
culate ujm over the interval < M < 5. This covers all sub- 
sonic (hot) and a wide range of supersonic (cold) discs. We 
calculate the normal modes for a=0.5, 0.1 and 0.03, which 
respectively require J = 50, 100 and 200 for a computational 
accuracy of 1%. 

Top row in Fig. 2 shows the eigenfrequency spectrum 
of the m = mode for M=0.7, 2.5, 5 and a = 0.5. Most 
eigenfrequencies of the model with M=0.7 are real (Fig. 2a). 
There are a few divergent complex eigenfrequencies whose 
magnitudes are increased by increasing J. The complex 
eigenfrequencies completely disappear for 1.8 < M < 2.7 
(Fig. 26). A convergent complex branch occurs for M > 2.7 
whose eigenfrequencies correspond to non-rotating unstable 
ring modes with = (Fig. 2 c). The fundamental eigen- 
frequency, denoted by coc, is located at the end of the com- 
plex branch of the spectrum. Fig. 3 a shows the variation of 
the fundamental growth rate with respect to M and a. It is 
seen that for M > 3 our results do not vary that much as we 
change a. The upper limit of M in stable discs is Mcv = 2.7 
for a = 0.5 and Mcr = 2.43 for a = 0.03. Our results for Mcr 
can be compared with the critical v/a computed from equa- 
tion (6.2) of Lemos et al. (1991). Their v and a for the /3 = 
models are our 14 and co, respectively. By keeping the mi- 



nus sign before the square root in equation (6.2) of Lemos et 
al. (1991) (this is to guarantee ring formation), and setting 
L = 1, /3 = and 7 = 1, one finds v/a = 2.414. The dif- 
ference between this critical number and Ma = 2.43 of our 
Q — 0.03 model is only 0.66%, which validates the accuracy 
of our matrix formulation based on the expansions intro- 
duced in equations (25) and (26). We note that the choice 
of q; = 0.03 is quite enough for understanding the behav- 
ior in the scale- free limit. We have also found that unstable 
breathing modes do exist in subsonic discs with a = 0.5 if 
M < 0.56. For a = 0.03 we find M < 0.87, which coincides 
with the upper limit oiv/a calculated in Lemos et al. (1991) 
for the existence of breathing modes. 

The eigenfrequency spectra of the m = 1,2 modes 
have three branches in the complex {Qm, Sm)-plane. Sub- 
sonic discs have no convergent unstable (complex) branch. 
They are stable under m = 1, 2 excitations (Figs 2d and g). 
A continuum of convergent real eigenfrequencies exists on 
the ilm-axis, which survives for almost all Al . although it 
may be shrunk to a small interval of positive pattern speeds. 
Aoki et al. (1979) have encountered similar eigenfrequen- 
cies in studying the gaseous Kuzmin disc. Convergent real 
eigenfrequencies correspond to pure oscillatory waves. Here 
we are interested in growing modes and we will overlook the 
real branch of the spectrum. As we increase M, a convergent 
unstable branch, with the largest arg(a)m)=arctan(sm/f^m), 
occurs in the spectrum. One can readily identify the funda- 
mental eigenfrequency Uc on this newly born branch. In a 
sufficiently cold disc, the convergent branch attracts more 
eigenfrequencies of the linear system (35) and Im(a'c) be- 
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comes the global maximum (see Figs 2/ and 2i). In Figs 3b 
and c we have illustrated arg(tJm)=arctan(sm/fim) {rn — 
1, 2) in terms of M and for several choices of a. These fig- 
ures show that the growth rate tends to zero faster than 
the pattern speed while the disc is stabilized. Instabilities 
are suppressed in discs with M < 1 because pressure waves 
travel faster than gravitationally generated waves. Figs 36 
and c can be compared with Fig. 3 in GE. The graphs of 
aig{LOm) have similar trends in our and their papers: For 
m = 1, arg(cJm) has a local minimum, while it is an almost 
monotonic function of M for cold discs and for m = 2. The 
major discrepancy between our and GE results occurs when 
M ^ 1. GE predicted a rapid growth of arg(ti;„i) contrary 
to the abrupt drop in our calculations. We ascribe the ex- 
isting discrepancies to the type of imposed internal bound- 
ary conditions. We use a cutout, which disables the medium 
needed for the propagation of incoming waves. GE, however, 
fix arg(tj) to make a reflective boundary at i? = 0. 

Fig. 4 displays the unstable rn = 0, 1 and 2 modes of 
a circular disc with M = 5 and a = 0.5. We have plotted 
the positive part of Am,o{R) cos[m(j) + Qm,o{R)]- Gray scale 
contours mark 10 to 100 percent of the maximum perturbed 
density. The m — 1 and m — 2 modes are tightly wound spi- 
rals whose amplitude functions, Aifi{R) and A2fi{R), have 
only one peak. There is no local density concentration along 
the spiral arms and the wave amplitude slowly decays well 
within the corotation radius. We have also plotted ver- 
sus j below each mode shape in Fig. 4. The plots show a 
rapid convergence of the sequence |a™|. As it can be seen 
in mode shapes, unstable waves cannot penetrate into the 
centre. This is due to the inner cutout, which damps the im- 
pinging waves by immobilizing the mass of central regions. 



5 MODES OF NON-AXISYMMETRIC DISCS 

We now apply the formulation of §3 and calculate the m = 0, 
1 and 2 modes of our non-axisymmetric discs. We set L = 3, 
which keeps nL + 1=7 interacting modes. Such a trunca- 
tion has enough accuracy for the models with M > 1, which 
mostly have e < 0.24 (see Fig. 1). The major limiting factor 
in the normal mode calculation of non-cixisymmetric discs 
is the choice of a whose small values need too many terms 
in (25) and (26) to ensure the convergency. Taking into ac- 
count the fact that the size of A in non-axisymmetric discs 
increases by a factor of (nL -|- 1) x (nL -|- 1) compared to the 
axisymmetric limit, the run time of our eigenvalue and eigen- 
vector solver rises drastically. To optimize the cost of numer- 
ical computations, we use a — 0.5 in our non-axisymmetric 
discs because it requires only J = 50 to achieve a relative 
accuracy of 1%. This is a legitimate choice because accord- 
ing to the graphs of Fig. 3 the qualitative features of normal 
modes do not change by varying a. 

Fig. 5 shows the evolution of the fundamental eigonfro- 
quencies with respect to the variations in e for two choices 
of M for each mode number m. Left panels in Fig. 5 corre- 
spond to a marginal M, which is chosen slightly larger than 
the upper limit of Mach number for stable circular discs. 
Right panels in Fig. 5 arc for a cold disc with M = 5. Left 
panels in Fig. 6 display the m = 0, 1 and 2 modes of non- 
axisymmetric discs that bifurcate from the modes of Fig. 4. 
Gray scale contours show the positive part of 



Ei'"'(7?,0,e) = A„(J?,0,e)cos[m0-he™(J?,0,e)]. (38) 

Right panel to each mode shape shows the absolute values of 
the scaled expansion coefficients, la*' |=e''' | {k = m + nl), 
in terms of I and j. 

Our computations show that non-axisymmetric discs 
with M < 2.7 are stable for m=0 excitations independent of 
the magnitude of e. In discs with M > 2.7 a tangent bifur- 
cation occurs in the complex {ilm, Sm)-plane at the location 
of each eigenfrequency (including lOc) and two sets of eigen- 
frequencies with non-zero pattern speeds come to existence. 
They are of the form a;o=±f2o + iso- This is an interesting 
result: we have found an unstable, ring-shaped density per- 
turbation that can rotate. Our computations show that the 
models with 2.7 < Af < 2.91 are stabilized by increasing 
e. Fig. 5a shows that how a model with M=2.8 is even- 
tually stabilized for e « 0.042. The pattern speed of the 
m = mode takes both positive and negative values because 
equation (23) is invariant under the reflections t ^ —t and 
(p —> —(j) for m = 0. According to Fig. 1, the upper limit of e 
decreases for our cold bisymmetric discs, and therefore, mod- 
els with large M cannot be stabilized only by varying the 
non-axisymmotry parameter. Nonetheless, the placement of 
eigenfrequencies in the complex plane is similar for all un- 
stable discs: the pattern rotates faster and grows slower as e 
increases. Fig. 6a displays the growing m=0 mode of a non- 
axisymmetric disc with (M, e) = (5, 0.015). The correspond- 
ing |a™+"'| (Fig. 6b) tend to zero as \l\ and j arc increased. 
Thus, the perturbation expansion given in (22) is convergent 
both in the azimuthal and in the radial directions. 

Our eigenfrequency calculations for m = 1 has been 
summarized in Figs 5c and d for M =1.7 and 5, respectively. 
The unstable m = 1 mode has a definite sense of rotation 
and its pattern speed decreases as e is increased. For the 
disc with M = 1.7, ujm steeply falls off and ultimately van- 
ishes for e « 0.05. This property is shared by all models in 
the range 1.6 < M < 2.3. Models with M < 1.6 are stable 
for m = 1 excitations regardless of the magnitude of e. The 
fundamental eigenfrequency of the models with M > 2.3 
does not vary considerably as e is increased in the admissible 
zone of Fig. 1. However, the mode shape is highly influenced 
by the non-axisymmetry. Fig. 6c shows the unstable m = 1 
mode of a non-axisymmetric disc with (M, e)=(5, 0.015). Al- 
though the overall shape is a single axmed spiral, which cor- 
responds to the I = wave component of equation (14), a 
number of dense clumps occur along the major spiral arm. 
These local clumps are generated by the |^| > wave compo- 
nents of the perturbed density expansion. Fig. 6d, shows the 
variation of |a™^"'| versus I and j for rn = 1. The sequence 
nicely converges as |/| is increased. It is seen that the choice 
of L = 2 provides enough accuracy both for m = and for 
m = 1. 

Finding the eigenfrequencies of the m = 2 mode for 
AI < 2 is troublesome because the convergent complex 
branch emerges very close to divergent branches (sec Fig. 
2h). We have written a program that computes the full eigen- 
value spectrum for several choices of J and searches between 
all eigenvalues for the most convergent one. According to 
our program, convergent unstable rn = 2 modes exist for 
M > 1.4. Numerical results displayed in Fig. 5e and / show 
that the fundamental eigenfrequency is robust to variations 
in e. The corresponding mode shapes are double armed spi- 
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il2 £lz 

Figure 5. Evolution of the cigcnfrcqucncics of non-axisymmctric discs witli respect to e and M. Top, middle and bottom panels 
correspond to m = 0, 1 and 2, respectively. Right panels correspond to M = 5. (a) M=2.8 (c) M=1.7 (e) M = 1.4. As Panels a and 
c show, the m = and m = I modes arc stabilized by increasing e in the allowable zone of Fig. 1. The m = mode can be stabilized 
for 2.7 ^ M < 2.91. This interval becomes 1.6 ^ M < 2.3 for m = 1. The eigenfrequencies do not vary substantially for cold discs with 
M = 5 and for the m = 2 modes. 



rals with some sort of dense clumps distributed along the 
major arms. The pitch angle of the global arms deceases as 
the disc is cooled (large M). Fig. 6e displays the m = 2 
mode of a non-axisymmetric disc with (M, e)=(5, 0.01). The 
perturbation series (22) slowly converges in the ^-direction 
for m = 2. So we kept more azimuthal wave components 
and used L = 3 (Fig. 6f). 



6 DISCUSSIONS 

The classical procedure in the instability analysis of gaseous 
(and stellar) discs is to start with an axisymmetric equi- 
librium and then perturb the governing equations in the 
vicinity of this equilibrium state. The subsequent eigenvalue 
problem [e.g., equation (35)] results in a complex frequency 
ujm = +iSm, which implies an m-fold barred/spiral pat- 
tern Em = Am.(R) cos[in(f> + Qm{R)] that rotates by the 
angular velocity Qm/m and exponentially grows by the rate 
Sm- Nevertheless, most events of physical significance hap- 
pen in nonlinear regime when the modes begin to interact. In 
such a circumstance, more than one azimuthal wave compo- 
nent is present in the response spectrum. The most general 
form of a density wave in nonlinear regime is 

-Hoo 

Si(Ji,0,t)= ^m(i?,t)e'["'^+®'"(^'*'l, (39) 

m= — oo 



which means that the amplitude and phase angle of Ei arc 
functions of all spatial and time variables. A primary result 
of (39) is the generation of dense clumps at the locations 
of the maxima of Ei. Clumps may undergo gravitational 
collapse, merge and give birth to gas giants. This is a widely 
accepted formation scenario of protoplanets (e.g., Boss 1997; 
Mayer et al. 2004) because it can take place in short time 
scales. Asymmetries like the presence of a double star (Boss 
2006) can also generate dense clumps. 

According to (31) and calculations of §5, unstable den- 
sity waves of non-axisymmetric discs have a rich structure 
(in linear regime) even for very small non-axisymmetry pa- 
rameter e. Although the temporal evolution of the mth mode 
is governed by the simple function exp(iciJmi), the eigenfunc- 
tion associated with uim is by no means simple and it con- 
tains infinite number of azimuthal wave components, which 
in turn are capable of generating local clumps. Since realistic 
accretion discs have not initially the nice axisymmetry fea- 
ture, their instability can give birth to protoplanets sooner 
than what nonlinear hydrodynamical simulations predict. A 
protoplanet, as a possible product of the gravitational col- 
lapse of a dense clump, is separated from the rest of the 
evolving gas and migrates until it is captured by a resonant 
island in the phase space. 

Our perturbation theory can be applied to any non- 
axisymmetric configuration regardless of the nature of the 
centre. Cutouts, which are special treatments designed for 
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Figure 6. Mode shapes (left panels) of a set of non-axisymmetric discs with M = 5.0 and a = 0.5. Top, middle and bottom rows 
correspond to (m, e) = (0,0.015), (1,0.015) and (2,0.01), respectively. These modes bifurcate from the modes of Fig. 4. Gray scale 
contours mark 10 to 100 percent of the maximum perturbed density with a step of 10. Local density peaks (darker regions) along the 
major spiral arms correspond to the local maxima of Sj™' (K, e). Right panel of each mode shape shows the corresponding absolute 
value of aj = el'la* (k = m + nl) versus j and I. The results clearly show the convergence of our perturbation series both in the radial 
and in the azimuthal directions. 
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scale- free models, are not needed in cored systems with finite 
total mass. The most important factor in success of our per- 
turbation theory is the smallness of the non-axisymmetry 
parameter e, which prohibits the interaction of azimuthal 
wave components in the zoroth order terms. Our equilib- 
rium density has a single azimuthal Fourier component [see 
equation (3)]. Accordingly, the harmonic numbers of wave 
components differ in steps of nl {I = 0, ±1, ±2, • • •). For 
bisymmetric discs studied in this paper, this means that 
wave components have only oven or odd harmonic numbers. 
We will have the most complete spectrum if the equilibrium 
fields include Fourier terms with both even and odd har- 
monic numbers, or if the disc is lopsided with n — 1. 

We presented three length scales to the perturbed equa- 
tion (14). The first two are the cutout radii, 7?in and Rout, 
and the third one is the length scale of the basis func- 
tions, Rc- At a first look, such a treatment of scale- free 
systems seems unfavorable because in this way we destroy 
their scale-invariance and cuspy features. However, with ap- 
propriate selection of the exponents of the cutout functions, 
A'in=A'out=2, and the weight function of the velocity com- 
ponents, W2(C)) the coefficient matrix A of our eigenvalue 
problem became a function of the two dimensionless param- 
eters a and 7. We easily decreased 7 and noticed that un- 
stable waves do not extend to distant regions although the 
mass of our equilibrium models (with flat rotation curves) 
is infinite at 7? ^ 00. We then let a tend to zero with the 
aim of removing the inner cutout. This needed more terms 
of the series of associated Legendre functions, but the qual- 
itative features of our results remained unchanged as long 
as mode shapes and arg(ciJm) were concerned. There was an 
excellent match between our results and those of Lemos et 
al. (1991) who used logarithmic spirals to model their per- 
turbed quantities. Logarithmic spirals do not introduce any 
length scale to variational equations. 

The critical value of M, below which our bisymmet- 
ric discs are stable, is a function of the fundamental wave 
number m. Numerical results show that the models with 
a = 0.5 are stable for m = excitations when M < 2.7. 
This limit approximately becomes 1.6 and 1.4 for m = 1 
and 2, respectively. The most interesting property that we 
have found for unstable non-cixisymmetric discs is that the 
m = and m = 1 modes can be stabilized by increasing 
the non-ajcisymmetry parameter e. This happens for m = 
when 2.7 < M < 2.91 and for m = 1 when 1.6 < M < 2.3. 
There are some remarkable differences between the stabi- 
lized m = and m = 1 modes. The stabilized m = mode 
has the form of a ring, which rotates according to the results 
displayed in Fig. 5 a. It occurs when the wave components 
with the azimuthal harmonic numbers m — nl — —nl and 
m+nl = nl cancel each other for I > 0. The stabilized m = 1 
mode, however, is a non-spiral, equiangular neutral mode. It 
emerges for a sufficiently large e, when the wave component 
with the harmonic number m — n = —1 {I = —1) cancels 
the fundamental component with m -I- n x = 1 (/ = 0). It 
is worth noting that the m = 2 mode cannot be stabilized 
by increasing e, because the interacting wave components of 
0{e) have the harmonic numbers m — n = {I = —1) and 
m -I- n = 4 (/ = 1). Neither of these first-order components 
is capable of suppressing the effect of the fundamental com- 
ponent with the harmonic number m + n x = 2. On the 
other hand, the wave component with the harmonic number 



m — n X 2 = —2 is of 0{e^), which is too weak to com- 
pensate the zeroth order m = 2 component. Based on the 
above reasonings, we speculate that the m = 2 mode can be 
stabilized in discs whose surface density includes the fourth 
harmonic of the azimuthal angle or when n = 4 in the simple 
discs of JA. Furthermore, the unstable m = 1 and m = 2 
modes of lopsided discs (with n = 1) cannot be stabilized 
by increasing the non-axisymmetry parameter. 
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APPENDIX A: THE ELEMENTS OF THE 
LINEAR OPERATOR L 

The linear operator L appeared in equation (17) is in terms 
of the active surface density Eact(^, 0), and the velocity field 

{U20,U3o)^[vRo{R,<j)),V^o{R,4>)] , (Al) 

of the equilibrium state. The elements of L are 

ill = =; Ot + + VroOr + -^04, 

^act \ it it 



L12 = ;^ + [afl-f-afl(lnE,ct)], 



(A2) 
(A3) 
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Ll3 = -^[dcl> + 84, (In Sact)] , 



[ai{-9ij(lnEact)], 



L22 = dt + dnVRo + vroOr + ^^84,, 
1 2u^() 

-L23 = -^O4,VR0 

-£-32 = dRV4,o + ^ , 

^33 = dt + -^d^v^o + ^ + vroBr + ^9^,. 



(A4) 

(A5) 
(A6) 
(A7) 

, (A8) 
(A9) 

(AlO) 



APPENDIX B: THE ELEMENTS OF THE 
OPERATOR MATRICES Dp 

The operators 'Dp=[Dpaj] = 1,2,3) appear in the in- 
teracting terms of (23). They arc functions of R and d/dR. 
Let us define k = m + nl and introduce 



H{R) 



R^ 



H,(«) = .+oV+3.»(|-) -<.'(f) , 

' 5M^ - 1 + 3n{M^ - 1)' 



/Ui = 1 - 



M2 



2(1 + n) 

4 + Se^ + 8n(l - 2n) (l + e^) 
" 8(l + n)(-l + 4n2) 
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2 + 71 ^ 72n^ + 104n2 
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- 1 
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4(l + n)(-l + 4n2) 

.2 12n^ + 18n2 - 5n - 5 
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-16n^-8n^ + l 
8(l + n)(-l + 4n2) 



+M- 



2 24n^ + 32n^ - lOn - 9 



8(l+n)(-l + 4n2) 
For IpI < 3, we have found 
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